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Abstract. This work is on the numerical approximation of incoming solutions to Maxwell's 
equations with dissipative boundary conditions, whose energy decays exponentially with time. Such 
solutions are called asymptotically disappearing (ADS) and they play an important role in inverse 
back-scattering problems. The existence of ADS is a difficult mathematical problem. For the exterior 
of a sphere, such solutions have been constructed analytically by Colombini, Petkov, and Rauch |10| 
by specifying appropriate initial conditions. However, for general domains of practical interest (such 
as Lipschitz polyhedra), the existence of such solutions is not evident. 

This paper considers a finite-element approximation of Maxwell's equations in the exterior of 
a polyhedron, whose boundary approximates the sphere. Standard Nedelec-Raviart— Thomas el- 
ements are used with a Crank-Nicholson scheme to approximate the electric and magnetic fields. 
Discrete initial conditions interpolating the ones chosen in |10| are modified so that they are (weakly) 
divergence-free. We prove that with such initial conditions, the approximation to the electric field is 
weakly divergence-free for all time. Finally, we show numerically that the finite-element approxima- 
tions of the ADS also approximates this exponential decay (quadratically) with time when the mesh 
size and the time step become small. 
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1. Introduction. This paper studies the numerical approximation of incoming 
solutions to Maxwell's equations in terms of the electric field, E(t,x), and the mag- 
netic field, B(t, x), 

eE t - curl yT x B = -j, B t + curl E = 0, 
dive.E = 0, divB = 0, 

in the exterior of a spherical obstacle, with dissipative boundary condition on the 



sphere (see (1.1)). Here, e is the permittivity of the medium, fj, is the permeability, 
and div j — 0, where j is the known current density of the system. We approximate 
numerical solutions, whose total energy decays exponentially with time. Such solu- 
tions are called asymptotically disappearing (ADS) and this phenomenon is of interest 
for inverse back-scattering problems, since the leading term of the back-scattering 
matrix becomes negligible. Details and construction of such solutions for the exterior 
of the unit sphere are found in the recent work by Colombini, Petkov, and Rauch [TP] . 
The dissipative boundary conditions of interest are 

(1.1) (1 + j)E tan = —n A fi~ 1 B tan on the boundary of the obstacle. 

Here, n is the outward unit normal to the boundary and 7 is a parameter satisfying 
7 > 0. 
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For a spherical obstacle, a case on which we focus on here, the boundary is \x\ = 1. 
In [TU] , it is shown that for any value of the parameter 7 > determining the dissipa- 
tive boundary condition, there exist initial conditions such that the boundary value 
problem for Maxwell's equations has a solution, which decays exponentially in time 
as 0(e rt ), with r < 0. It is also interesting to note that in space such solutions also 
decay asymptotically at infinity, i.e. they behave as 0(e r ^) (see [2D]). Moreover, for 
dissipative boundary conditions, (1.1 1, if 7 > 0, there are no disappearing solutions, 
u(T, x), that vanish for alH > T > in the exterior of the sphere (see |12j). 

The focus of this work is on the finite-element approximation of the ADS in the 
exterior of a polyhedron that approximates the sphere. This is a first step towards 
developing numerical techniques, which later can be used to construct approxima- 
tions to ADS for more complicated obstacles and also for other symmetric hyperbolic 
systems with dissipative boundary conditions. Such a general study is feasible due 
to a recent result of Colombini, Petkov, and Rauch [11] for hyperbolic systems whose 
solutions are described by a contraction semigroup V(t) = e Gt , t > 0. More precisely, 
it was shown in [llj that if certain coercivity estimates are satisfied, then, the spec- 
trum of the generator, G, in the left half plane, Re(z) < 0, is formed only by discrete 
eigenvalues with finite multiplicities. Every such eigenvalue, A : Re(X) < 0, yields an 
ADS solution, u(t, x) = e At /(x), with Gf = Xf. On the other hand, the existence and 
the location of such eigenvalues for less regular obstacles is a difficult and interesting 
mathematical problem (both from theoretical and numerical point of view) , albeit it 
falls beyond the scope of the work reported here. 

For the finite-element spaces that we use, their properties and implementation in 
the numerical models based on Maxwell's system are more or less known. Classical 
references on the piecewise polynomial spaces relevant in such approximations are the 
papers by Raviart and Thomas [21] (for two spatial dimensions) , Nedelec [TU [17] , and 
Bossavit 0. The method that is used here is an application of the techniques devel- 
oped by Brezzi [8., (see also Brezzi and Fortin [5]). Many results and references on 
Maxwell's system and its numerical approximation are found in Hiptmair's work [13] 
and in a monograph by Monk [T5]. In most of these works, the emphasis is on systems 
with perfect conductor boundary conditions. This important difference between our 
work and the previous ones is that, here, we apply the methods to a problem with 
dissipative boundary conditions, where the electric field, E, and the magnetic field, 
B, are paired together. A related discussion on the role of the boundary conditions in 
the numerical solution of systems of PDEs are discussed in a recent paper by Arnold, 
Falk and Gopalakrishnan [2J. 

We consider the finite-element problem associated with Maxwell's equations in 
a finite domain between two spheres, £1 = {x | 1 < \x\ < R}, for a fixed, large 
enough R. The discretization of the Maxwell's equations that keeps the coupling 
between the electric and the magnetic field in tact and can be derived via the mod- 
ern techniques in finite element exterior calculus is described in Arnold, Falk and 
Winther [TJ [3] . We approximate the region by first constructing a tetrahedral mesh 
for the domain = [-R/2, i?/2) 3 \[— 1/2, 1/2] 3 and then mapping it to £1 in polar 
coordinates, Qx\t 3 , 6 , (j>) H> (1x^,0,0). We denote the resulting partition in tetra- 
hedra by Tk- I n such a way, we obtain a computational domain (denoted again with 
f2), which is polyhedron and which approximates the exterior of the annular domain 
{x I 1 < I a; I < R}. Note that the ADS constructed in [TU] for the exterior of the 
sphere do not need to satisfy the dissipative boundary condition on the polyhedron. 
However, we show numerically that the finite-element solution with initial conditions 



Numerical Approximation of Asymptotically Disappearing Solutions of Maxwells Equations 3 



approximating those in |10j yields good approximation of the ADS constructed an- 
alytically, when the mesh size becomes small and the polyhedron gets closer to the 
sphere. 

This paper is organized as follows. In Section [2j we introduce some notation and 
state the strong form of the boundary value problem for Maxwell's equations that 
we consider. Section [3] describes the variational (weak) formulation and discusses 
the energy decay of the corresponding system. Next, we discuss the discretization of 
this variational form and how we can guarantee a good approximation of the ADS in 
Section[4j In Section[5j we describe the matrix representation of the semi-discrete sys- 
tem and the properties of the Crank-Nicholson scheme that we use for time stepping. 
Numerical results for the sphere are shown in Section [6j Finally, concluding remarks 
and discussions on constructing initial conditions as well as a choice of parameters for 
more complicated obstacles are presented in Section [7j 

2. Notation and preliminaries. First, some standard notation is introduced, 
which is needed in the following sections. The Euclidean scalar product between two 
vectors in IR rf is denoted by (a,b) and the corresponding norm is \a\ 2 = (a, a). The 
standard L 2 (f2) scalar product and norm are denoted by (•, •) and || • ||, respectively. 

2.1. Maxwell's system. The system of partial differential equations (PDEs) 
of interest is Maxwell's system with a dissipativc boundary condition (impedance 
boundary condition). Let O be a bounded, connected (could be convex) domain, 
O C 1R 3 . As stated in the introduction, we consider Maxwell's equations in the 
exterior of O, that is, in Q = R 3 \ O, which is as follows: 

(2.1) eE t -cmln- 1 B = -j, 

(2.2) B t + curl E = 0, 

(2.3) dives = 0, 

(2.4) divB = 0. 

For the rest of the paper, we assume that the permitivity, e, and the permeability, 
/i, are equal to 1. Future work will involve investigating numerical methods when these 
parameters are allowed to vary with the domain. We also assume that f2 is bounded, 
which means fl = S\0, where S is a ball in IR 3 with sufficiently large radius. While 
in general this could be a restriction, in this case it is not, since the solutions that are 
approximated decay exponentially when |x| — > oo. 

2.2. Dissipative boundary conditions. The boundary conditions that are of 
interest are also known as impedance boundary conditions. To state such type of 
boundary conditions, we first define 

Ti = dttndO, r = dn\r l , 

where Ti represents the boundary of the inner obstacle, and r o the outer boundary 
(i.e., the boundary of S). The orthogonal projection, Qtam on the component of a 
vector field tangential to I\ is also needed, which for any x E Ti and a vector field 
F(x) g IR 3 is defined as the tangential component of F(x), namely: 

F tan = Qtan-F = F - (F, n)n = -n A (n A F) , 

where n is the normal vector to the surface Ti . We denote by n the normal pointing 
away from the domain (i.e. pointing into O). We note that all the quantities above 
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depend on s e Tj. The boundary condition of interest is the one in (1.1) and it is 
recalled here: 

(1 + 7)-E tal i = —n A B taa or equivalently (1 + 7)-E tan = -n A B. 

As pointed out above, 7 > is a constant, i.e. 7 € H. However, the same methods 
can be applied for 7(2;) > as a function on the boundary of the domain. 

Remark 2.1. Note that for a perfectly conducting obstacle, the tangential com- 
ponent of E vanishes on the boundary, namely: 

E An = 0, xedtt. 

However, again, this paper considers the case of an impedance condition, where the 
obstacle is not a perfect conductor. This is closer to real-world applications, where 
dissipative boundary conditions occur frequently. 

3. Function spaces and variational formulation. 



3.1. Function spaces. To approximate the differential problem (2.1 )-( 2.4 ) with 



the boundary conditions given in (1.1), the function spaces for the problem at hand 
need to be identified. Given a Lipschitz domain fl and a differential operator T>, a 
standard notation for the following spaces is used: 

H(p) = {ve (L 2 (n)) d ,Vv e L 2 (n)}, 

with the associated graph norm, 

\\uf V;n =\\uf + \\Vuf. 

By taking T> = div or T> = curl, the Sobolev spaces TJ(div) and H(cur\) are obtained. 
Also, notice that 

H\n) = ff(grad) , L 2 (il) = H(id) . 

For example, if (curl) is the space of L 2 (f2) vector- valued functions, whose curl is also 
in L 2 (f2). Similarly for iJ(grad) and if (div). 

The following three spaces are needed (the first one for scalar functions and the 
second and third for vector- valued functions): 

#o(grad) = = {v € H\n) such that v\ gn = 0}, 

-ff imp (curl) = {v e ff(curl) such that v A n\ =0}, 
H i m p (div) = {»£ H(div) such that (v, w) L =0}. 

Note that the tangential component of the electric field E and the normal component 
of the magnetic field B on the outer boundary, T D — dil \ T i} are set to zero for the 
elements of ifi mp (curl) and i/i mp (div). 

Another issue to address is related to the fact that the boundary of the compu- 
tational domain consists of two connected components F\ and r o . In such a case, the 
electric field E is unique up to a harmonic form (harmonic function, constant on Ti). 
To resolve the ambiguity, we consider electric fields in a subspace of H lrnp (curl) , which 
is orthogonal to the one-dimensional space of harmonic forms. Let f) be the unique 
solution to the Laplace equation: 



(3.1) 



-Afj=Q, f) = l on Ti, and f) = on F f 
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Then, define -ff; mp (curl) as the space of functions orthogonal to grad fj: 

(3.2) i?im P (curl) = {v £ -ffimp (curl) such that (v, gradf)) = 0}. 

Finally, for the time-dependent problem considered here, the relevant function spaces 
are 

iJ (grad;£) = {v(t, •) G H^(Q) for all t > 0}, 
fli mp (curl; t) = {v(t, •) £ iJ imp (curl), for all t > 0}, 
-ffim P (div; t) = {v(t, •) e flimp(div), for all t > 0}. 

In another words, if H(T>) denotes any of the Hilbert spaces H (gra,d), ff; mp (curl), 
or i?; mp (div), then H(T);t) is the space of functions, which for each t £ [0,oo) takes 
on values in H(T>). We assume that the elements of any of the spaces Ho(gi&d;t), 
(resp. -ffi mp (curl; t), or _ffi mp (div; t)) are diffcrcntiablc with respect to t as many times 
as needed. We refer to McLean [14] and also Monk [15] for properties of the above 
spaces, and related density results. 

3.2. Variational formulation. Next, we derive a weak form which was shown 
to us by D. N. Arnold [I]. We introduce p € iJo(grad;i) such that 

(3.3) (pt, q) = {E, gradg), for all q e i? (grad), p(0, x) = 0. 

This is an auxiliary variable associated with the constraint that E is divergence-free. 
If the initial condition for E is divergence-free, then as shown below, p is zero for all 
times. However, if the initial guess is not divergence-free, then p does not have to be 
zero. 



First, multiply equations (2.2 1 and (2.1| by test functions C £ ifj mp (div) and F e 
^imp(curl), respectively. Next, integrate by parts and use boundary conditions (1.1 1 
and the identities 



(n A B, F) = ((n A B tan ), F tan ), (E tan , F tan ) = (n A E, n A F). 
to obtain for all C £ i?i mp (div), and for all F £ i7j mp (curl), 

(B t ,C) + (cmlE,C) = 0, 

(3.4) (E t ,F) + (Vp,F)-(B,cm\F) 

l(l| 7 )/ ri (nA£, H AF) d 7 = -(j,F), 

Finally, we get the following variational problem: 

Find (E, B,p) £ 77; mp (curl; t) x _ff imp (div; t) x i/ (grad; t), such that for all (F, C, q) £ 
-ffimp(curl) x _ff imp (div) x Hq(£1) and for all t > 0, 

(3.5) (E u F) = - (gradp, F) + (B, curl F) - + (E tan , F tan ) - (j, F), 

(3.6) (B t ,C) = -(cwlE,C), 

(3.7) ( Pt ,q) = (E,gmdq). 

At t = 0, the following initial conditions are needed, 

(3.8) E(0, x) = Eq(x), B(Q, x) = B (x), p(0, x) = 0. 
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Here, E £ i? imp (curl), B £ if imp (div), and we assume that B is divergence-free. 
As it is well known, this assumption implies that B is divergence-free for all t > 0. 

Remark 3.1. Regarding the existence and uniqueness of the solution to the vari- 
ational problem, (3.5) -(3.7 1, we point out that in the case of a perfect conductor, 
solution results for Lipschitz polyhedra are found in Pauly and Rossi \18\/ and also 
in Birman and Solomyak [5, 6}. However, for the dissipative boundary conditions we 
consider and for Lipschitz polyhedral domains, to our best knowledge, there are no 
results on existence and uniqueness. For solution results on smooth domains (C ) we 
refer to the monograph by Petkov where the general case of symmetric hyper- 
bolic systems is considered and existence and uniqueness are derived using semigroup 
theory. For Lipschitz polyhedral domains it is not straightforward and difficult to see 
that the solution is related to a contraction semigroup. 

Next, the following proposition shows that for a divergence-free E £ iJ; mp (curl; t), 
the variational form satisfies the divergence-free condition for E weakly for all time. 

Proposition 3.1. Letu — (E,B,p) satisfy Equations ( 3.5 ) -( 3.7 ) and the initial 
conditions (3.8). If Eq £ Hi mp (cml) is weakly divergence-free, then, for all t and x, 
p(t, x) = and, hence, 



(E, grad q) = 0, for all q£H^(n). 



Proof. 

get 

(3.9) 



In Equation (3.7), take q £ Hq(Q), then differentiate with respect to t to 



(Ptt,q) = {E u gradg), for all q £ H^(Cl). 



In the first equation, (3.5), set F = gradg. Note that the tangential component of 



grad q is zero, because q £ H^ ip) has a vanishing trace on the boundary of il. Thus, 
the boundary term vanishes and we are left with 

(E t ,gradq) = - (grad p, grad q) + (B, curl grad q) - (j, grad q). 

Using integration by parts on the last term yields 

(j',gradg)= / (n • j)q - (divj, q) = 0, 

since j is divergence-free and q vanishes on the boundaries. Taking this into account, 



along with (3.9) and the identity curl grad = 0, rewrite Equation (3.5) as follows: 

(3.10) ( ftt , (? ) + (gradp,grad< Z ) = 0, for all qeH^Q). 

Since p(0) = and p t (0) = div E(Q) = 0, it is concluded that p(t,x) = for all t > 
and all x £ O, since it is a solution of the homogenous wave equation, i.e., Equation 
(3.10). The condition pt(0) = div E(0) comes from integrating (|3.7[) by parts and 



Then, from (3.7), the desired result for E 

□ 



noting that q vanishes on the boundary, 
follows. 

Next, we show an energy estimate. 

Proposition 3.2. Letu — (E,B,p) satisfy Equations (3.5) -(3.7) and the initial 
conditions (3.8). Assume that j = (i.e. 



no external forces). Then, the following 



estimate holds for all T > 0: 



>\p(T,-)\\ 2 + \\E(T,-)\\ 2 + \\B(T,-)\\ 2 <\\E \ 



\B4 2 
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Proof. Fix t and take (q, F, C) = (p, E, B). Summing up the three equations (3.5 1 



(3.7) gives 

(3.H) j t (Ml 2 + ll^ll 2 + ||B|| 2 ) = -2(1 + 7)11^11^(1.0 

This identity holds for any t and the proof is concluded after integrating with respect 
to time. □ 

4. Finite-element discretization. 



4.1. Domain partitioning. To devise a discretization of (3.5 H3. 7), we ap- 
proximate the exterior of the sphere by a polyhedral domain, which is decomposed as 
a union of simplices (tetrahedrons). The polyhedral domain and its splitting is ob- 
tained by mapping a corresponding splitting of a cube to a polyhedron with vertices 
on the sphere. 

We consider a cube f2 = (— R/2, R/2) 3 and we split it into -p- cubes each with side 
length h. Here h = 2~ J , for some J > 2. From this partition, we remove all cubes that 
have nonempty intersection with the open cube u> = (—1/2, 1/2) 3 . Finally, we split 
each of the cubes from the lattice into 6 tetrahedrons and this gives a partitioning of 
Q\ui into simplices. Note that this partition has the same vertices as the lattice and 



Sl\u = U KeTh K. 

From this, we obtain a polyhedron approximating the region between two spheres by 
mapping every vertex of the lattice with polar coordinates (\x\e 2 ,9, <j>) to flas^, 6, <j>). 

Clearly this maps the interior boundary of to the unit sphere, and the outer 
boundary to a sphere with radius R. An example is shown in Figures |4. 1^4.21 We 
note that when h —¥ the corresponding polyhedron approximates the region between 
the spheres. 




Fig. 4.1. Cube Obstacle Fig. 4.2. Sphere Obstacle 

4.2. Finite-element spaces. For the discrete problem, we use standard piece- 
wise linear continuous elements paired with Nedelec jTBl [T7] finite-element spaces. We 
introduce the spaces i?/ li i mp (curl) C i?i mp (curl), and Hh; lmp (div) C i?i mp (div) and the 
FE solution is denoted by (E h , B h ,p h ) e H htimp (curl)xH h ^ mp (div)xH h (gr&d). More 
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precisely, we define 

H hMp (div) = {C e H hMp (div), C\ K = a K + /3 K x, MK e %}, 
H hfi (gra,d) = {q e H^(tt), q\ K is linear in x, VK E %}■ 

The space i?/ lj i mp (curl) is a properly chosen subspace of Hh. imp (curl), which is orthog- 
onal to the gradients of the discrete harmonic forms (but not necessarily to gradf)). 
This intermediate space is defined as 

fffe, imp (curl) = {F e # imp (curl), F\ K = a K + b K A x 7 6 77,}. 

Next, the discrete harmonic form, f) , is defined as the unique piecewise linear con- 
tinuous function satisfying 

(grad t) h , grad q) = 0, for all q € Hft, (grad), 
f) /l = l on r,, 
t) h = on r . 



Then, 

^Mmp(curl) = {F e H h<imp (cuil), {F, gradf)' 1 ) - 0}. 

In the definitions above, a,K,bx € 1R 3 are constant vectors for every simplex K in 
the partition and (3k € H. 

Spaces corresponding to the time-dependent problem are analogously defined us- 
ing the definitions from the previous section. We denote these spaces by Hh i mp (curl; t), 
Hh,i m p(div;t) and Hh,o(grad;t), respectively. Note that v S i?i mp (curl; t) implies that 
the tangential components of v are continuous. The other spaces induce certain com- 
patibility conditions as well. For example, the requirement B € if (div; t) in the 



definition of the space 7J/ lj i mp (div; t) at the beginning of Section 4.2 is equivalent to 
saying that the normal components of the elements from .i mp (div; t) are continu- 
ous across element faces. It is also easy to check that q £ Hh,o(grad;t) C Hq(£1) in 
the definition above implies that q is a continuous function (because it is a piecewise 
polynomial function, which is in Hq(Q,)). 



Finally, it is important to note that Equations ( 3.5 H 3.7 1 make sense for all E 6 
H (curl). As stated above, the piecewise polynomial functions on tetrahedral partitions 
of f2 are in iJ(curl) if their tangential components across the faces are continuous. 
Such functions, however, do not necessarily have continuous normal component across 
the faces of the tetrahedrons. Thus, the approximation E h to E is not in if(div) even 
though E G ff(div). 

4.3. Discrete weak form. After constructing the approximating spaces, the 
discrete problem is constructed by restricting the bilinear form onto the piecewise 
polynomial spaces. In the following, we set j — because we are interested only in 
the dependence of the solution on initial conditions. Denoting 

H h = i?/ l ,i mp (curl) x H htimp (div) x H hi0 (grad), 
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and restricting (3.5)-(3.7) to Hh leads to the following approximate variational prob- 
lem: Find (E h , iFV) € H h such that for all (F h , C h , q h ) G H h 

(4.1) (E' t i ,F h ) = -(gmdp h ,F h ) + (B h ,cuY\F h )-(l+-f) J (n A E h , n A F h ) 

(4.2) (B?,C h ) = -(cm\E h ,C), 

(4.3) {p h t ,q h ) = {E\g^Aq h ). 

In the following, the superscript h is omitted, since the considerations in the rest of 
the paper are focused on the discrete problem in H^. We also interchangeably use u 
and (E,B,p) and, similarly, w and (F,C,q). 



An operator is introduced, such that A : Hh *-> Hh via the bilinear forms in (4.1 1 



( |43| ). For u = (E,B,p) e H h and w = (F,C,q) e H h , set 

(Au,w) = -(E,gradq) + (gradp,F) - (B,curlF) + (cm\E,C). 

Corresponding to the boundary term, we also have the operator associated with the 
impedance boundary condition, 



(Zu, to) = (1 + 7) J (n A E.n A F 



Since we are now on a finite-dimensional space, we write the semi-discrete problem 
(discretized in space and continuous in time) as a constant coefficient linear system 
of ODEs, i.e. 

(4.4) u=-(A + Z)u 

From the definitions of A and Z it is obvious that A is skew symmetric and Z is 
symmetric and positive semi-definite. 

5. Matrix representation and time discretization. We now show that the 



assembly of system (4.4) can be constructed using only mass (Gramm) matrices 
formed with the bases in Hh, imp (curl), £f/ii mp (div), and Hh,o(grad) spaces and in- 
cidence matrices, whose entries encode the relationships "vertex incident to an edge", 
"edge incident to a face", etc. The aim of this section is to provide some insight 
into the implementation of such finite-element schemes and also to set the stage for 
presenting the Crank-Nicholson discretization in time. 

5.1. Matrix representation. We start with a description of the standard (canon- 
ical) bases in Hhfi(grad), -ff/i,i m p(div), and Hh, i mp (curl), respectively. By boundary 
vertices, edges, and faces we mean vertices, edges, and faces lying on the boundary, 
dft. For an edge, this means that both its end vertices are on the boundary and for 
a face it means that all three of its vertices are on the boundary. The remaining ver- 
tices, (edges, faces) are designated as interior vertices (edges, faces). We note that by 
a standard convention, it is assumed that for the triangulation in hand the directions 
of vectors tangential to edges and normal to faces are fixed once and for all. It is easy 
and straightforward to check that a change in these directions does not change the 
considerations that follow. 

We then have the following sets of degrees of freedom (DoFs): 
• DoFs corresponding to the set of interior vertices {iCi}"=i : A functional (also 
denoted by Xi) is associated with an interior vertex X{ as Xi(q) = q(xt) for a 
sufficiently smooth function, q. 
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• DoFs corresponding to the set of all interior edges and all edges on IV For 
a sufficiently smooth vector-valued function, v, and an edge, e € £, the asso- 
ciated functional is e(v) = A J e v ■ r e , where T e is the unit vector tangential 
to the edge. The direction of the tangent vector, r e , is assumed to be fixed. 

• DoFs corresponding to the set of interior faces, T: For a sufficiently smooth 
vector- valued function, v, and a face / € J 7 , the associated functional is 
f(v) — tjt J ,v ■ rif, where rif is the unit vector normal to the face. 

As bases for the spaces H h ^(gra.d), iJ^, imp (curl), and H himp (div) we take the 
piecewise polynomial functions, which are dual to the functionals given above. For 
the space .ffj^ofgrad), we denote these functions by {fjYj^-i- They are piecewise 
linear, continuous, and satisfy X)~(<Pj) = 5k j, where 5k j is the Kroneker delta. 

The bases for the other two spaces iJ/ ti i mp (curl) and fffc,i mp (div) are then given 
in terms of the basis for ^^(grad). For an edge e € £ with vertices (xi,Xj) and a 
face / € T with vertices (xi, Xj, Xk), 

ip e = \e\(ipigrad(pj - (fj gradipt), 

£j = l/l(<Pi(g rad< /'j A grad^fc) - ipj(grad<p k A gradt^) + ip k (grad<Pi A grad 

Here, r e = (xj — Xi)/\xi — Xj\ and the ordering of (x^, Xj, Xk) in a positive direction 
is determined by the right-hand rule and the normal vector rif. We then have the 
following canonical representations of functions in i?/ li i mp (curl; t), -fffc, imip (div), and 
-ff^oferad): 



For functions that also depend on time, i.e. for the elements of ff? lj i mp (curl; t), 
-fffi,imp(div; t), and Hf l0 (grad; t), we have similar representations with coefficients de- 
pending on time as well. 

Remark 5.1. In the rest of the paper, the same notation is used for the functions 
from Hh and their vector representations in the bases given above. This is done in 
order to simplify the notation. 

The entries of the mass (Gramm) matrices for each of the piecewise polynomial 
spaces are then, 



Next, the following matrix representations of the operators defined in the previous 
section are introduced, 





(M e ) ee , = (lP e ,l/> e ,), (Mf) 



ff 



g, 



(grad<^,i/> e ), fC fe = (curl^, £/). 



The matrix form of (4.4) is now rewritten as follows. 



(5.1) 
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5.2. Time discretization. To discretize system (5.1| in time, a Crank-Nicholson 
scheme is used. We look at a time interval, t £ [0,T], and approximate the solution 
at t = kr, k = 1, . . . , — , with r a given time step. Let u k = (E k , Bk,Pk) be the 
discrete approximation at the current time t = kr, and u k -\ = (E k -i, B k -i,Pk-i) T 
be the approximation at the previous time t = (k — l)r. Then, the Crank-Nicholson 



formulation of (|5.1|) is 

'Ml 

M f 

M v/ 



-M{u k -u k -x) = -\{A + Z)){u k +Uk-x), whereX=l Mf 

T 2 



Rearranging the terms, we get the following linear system for the approximate solution 
at time step kr in terms of the solution at time (k — l)r. 

(5.2) (±M + + Z)\ u k = (jM + Z)) ufc-x. 

Next, we show that if the initial condition is weakly divergence-free, as in the 
continuous case, we have that pk will remain zero for all time and, thus, Ek is weakly 
divergence-free for all k. This is the discrete analogue of proposition |3.1[ 

Lemma 5.1. Assume that (E ,gr&dq) — for all q € (gr&d) . For the 



Crank- Nichols on scheme described in (5.2 1, pk = for all k and (Ek,gra,dq) = for 
all q £ -ff^oferad) and all k. 

Proof. Start with po — and Eg being weakly divergence-free. It is shown 
that: 

If p k = and G T M e E k = 0, then p k+ i = and Q T M e E k+ i = 0. 

This is the matrix representation of the assumptions and claims in the lemma. Setting 
a = 2/r and using the defining relations for the Crank-Nicholson time discretization, 



(5.1) and (5.2), the following linear system for E k +i, B k +i, and p k +i is obtained: 



aMeEk+i - JC T M fB k+x + M e Gp k +i + ZE k+ i 

(5.3) = aM e E k + K T MfB k — ZE k , 

(5.4) M fJCEk+i + aM f B k+1 = -M fKE k + aM f B k) 

(5.5) - g T M e E k+ i + aM vPk+1 = 0. 

Here, Z indicates the reduced matrix of Z applied to just the edge DoF, E. Multi- 



plying Equation (5.3) from the left by Q T yields, 

ag T M e E k+ i - Q T K T M f B k+1 + G T M e g P k+i + G T ZE k+ i 



= aG T M e E k + Q T K T M f B k - G T ZE k . 



Next, note that ICG — (or equivalently G T IC T = 0), since the curl of a gradient is 
zero. Also, since any q £ Hf l> o(gisud) C Hq(£1) is zero on Fj and r o and the tangential 
component of its gradient is zero on the boundary edges, then ZQ = (or equivalently 



Q Z — 0, since Z is symmetric). Thus, (5.3) simplifies to, 



ag T M e E k +i + g T M e g P k+i = o. 
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Adding this to a times Equation (5.5), then gives 



(5.6) (Q T M e g + a 2 M v ) Pk+1 =Q 

The above relation is the matrix representation of the variational problem: 
(gradp fe+ i,gradg) + a 2 (p k+1 ,q) = 0, for all q € iIfo,o (grad) , 



and taking q = pu+i then gives that Pk+i — 0. Finally, from this fact and using (5.5), 
it is immediately shown that Q T A4 e Ek+i = 0, concluding the proof. □ 
Thus, using the Crank-Nicholson scheme and appropriate initial conditions, one 
can guarantee that the discrete approximation to the electric field will be weakly 
divergence- free for all time. 

5.3. Solution of the discrete linear systems. To solve the system, we look 
at the matrix corresponding to ^M. + | (A + Z), which is on the left side of (5.2). 



We have to solve the system with this matrix at every time step. Using the incidence 



matrices as in (4.4), this operator is written as 



i i i / *M e -JC T M f M e G\ i 

M + -(A + Z) = -\ M f K 2 -M f )+-Z. 
2 2 \-G T M e 2 -M v ) 2 

Since, the mass matrices, AA e , -Mf, -M V) are all SPD and Z is symmetric positive 
semi-definite and only contributes to the edge-edge diagonal block of the system, the 
entire operator can be made symmetric by a simple permutation. Multiplying on the 

left by J = | —I will yield the operator 

/i i \ ( hMe -\K T M } +\M e G\ j 

which is now symmetric. Therefore, the final system to solve is 

(5.7) J (^-M + + Z)j u n = J Q-M z )) u n-u 

and a standard iterative solvers such as MINRES can be applied. This is used in the 
test problems below. 

6. Numerical Results. Here, we perform some numerical tests by solving sys- 



tem (5.7) using the Crank- Nicholson time discretization described in the previous 
section. To test for decay in the energy of the solution, we start with initial condi- 
tions and boundary conditions of the form described in [ID]. We take as the domain, 
the area between a polyhedral approximation of the sphere of radius 1 and a polyhe- 
dral approximation of a sphere of radius 4. The inner sphere represents the impedance 
boundary and the outer sphere is considered far enough away (and it is for the solu- 
tions we approximate) that a Dirichlet-like perfect conductor boundary condition on 
the outer sphere is used. In other words, we prescribe E An, B ■ n, and p = on the 
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outer sphere. The exact solution (taken from [TQl Theorem 3.2]) is given as follows: 
(6.1) E t = 



(6.2) B t = e 



r(\ X \+t) 




Different values of 7 yield different values of r in solutions (6.1)-(6.2). Following [TU] . 
we have that (E*, B*,0) solves Maxwell's system with an impedance boundary con- 
dition on r,- and 



r = 1/2 (l - y/1 + 4/7) 



For the tests below, we take 7 = 0.05 (r = —4). 

6.1. Approximation of the initial conditions. Since the solutions given 
above are not in the finite-dimensional spaces considered, we take an initial condition 
Eq, which is based on the piecewise polynomial interpolant of the exponentially- 
decaying solution given in equations (ROT) at t — 0. In other words, we set 



]e(E,{0,x))xP e (x). 



We further correct Eq to get an initial guess that is orthogonal to the gradients as well 
as the gradients of the discrete harmonic form, gradf)' 1 . This is done in a standard 
fashion by projecting out these gradients. First, we find s £ ifo,fc(grad), such that for 
all q £ f/o,/i(grad), we have 



(grad s, grad q) = (E , grad q), E = E ,h - grad s 



CEo^gradf^) 
|| grad t) h \\ 2 



gradf)' 1 . 



As a result, Eq is orthogonal to the gradients of functions in iJo,/i(grad) and also 
to the gradient of the discrete harmonic form. We note that this orthogonalization 
requires two solutions of Laplace equation. It is straightforward to see that if the 
initial guess satisfies such condition, then the solution E satisfies this condition for 
all times. Finally, Bo is computed as Bq — ~KEq, 

6.2. Numerical results. We test the approximation to the asymptotically dis- 
appearing solutions on a grid with 728 (h = 1/8), 4,886 {h = 1/16), 35,594 (h = 1/32), 
and 271,250 (h = 1/64) nodes on the domain. The computational domain is shown 
in Figure |4.2| We run a MINRES solver on the Crank-Nicholson system that is 
symmetrized, Equation ( |5.7[ ), for 20 time steps using a step size r = 0.1. 
The results arc shown in Tables 



6.1 



6.4 



where we display the 



|L 2 (n) norms. The total energy of this system is given by H-En^m) 
Tables |6.1||6.4 show that over time the L2 norms of the electric and magnetic 
fields decay. For each mesh size, it appears that the energy reaches some steady-state 
value, where it does not decay anymore. Figure [6~2] shows that this final energy value 
(at time step 20), ||£'||^ 2 ^ n ^ + ||-B||^ 2 m^, decreases as h 2 , when h —> 0. Thus, as the 
polyhedron domain more closely represents the spherical domain, as in |10j . the total 
energy should decay to zero as expected over time. 



\E\ 

+ 1 



L 2 (Q) 



and 
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Step 


\\E\\l 2 (si) 


\\B\\ L2 (n) 





0.906 


0.559 


1 


0.830 


0.334 


2 


0.723 


0.155 


3 


0.614 


0.116 


4 


0.526 


0.172 


5 


0.450 


0.242 


6 


0.380 


0.304 


7 


0.316 


0.349 


8 


0.262 


0.379 


9 


0.233 


0.389 


10 


0.234 


0.384 


11 


0.248 


0.373 


12 


0.266 


0.358 


13 


0.284 


0.339 


14 


0.292 


0.322 


15 


0.288 


0.312 


16 


0.276 


0.307 


17 


0.262 


0.303 


18 


0.251 


0.299 


19 


0.248 


0.293 


20 


0.246 


0.287 



Table 6.1 
Sphere Obstacle. 7 = 0.05. h = 1/8. 



Step 


\\ E 1 £ 2 (£2) 


\\ b \\l 2 (si) 





0.553 


0.426 


1 


0.451 


0.266 


2 


0.349 


0.153 


3 


0.267 


0.093 


4 


0.200 


0.086 


5 


0.147 


0.103 


6 


0.117 


0.103 


7 


0.107 


0.099 


8 


0.104 


0.097 


9 


0.095 


0.101 


10 


0.088 


0.103 


11 


0.085 


0.104 


12 


0.084 


0.104 


13 


0.085 


0.102 


14 


0.086 


0.100 


15 


0.086 


0.099 


16 


0.086 


0.099 


17 


0.085 


0.099 


18 


0.086 


0.097 


19 


0.085 


0.095 


20 


0.085 


0.095 




Table C 


.2 



Sphere Obstacle. 7 = 0.05. h = 1/16. 



Step 


\\E\\L 2( n) 


\\B\\ L2 (n) 





0.425 


0.405 


1 


0.308 


0.262 


2 


0.214 


0.163 


3 


0.148 


0.102 


4 


0.103 


0.061 


5 


0.074 


0.044 


6 


0.057 


0.039 


7 


0.046 


0.037 


8 


0.040 


0.035 


9 


0.035 


0.034 


10 


0.033 


0.034 


11 


0.032 


0.034 


12 


0.031 


0.035 


13 


0.031 


0.034 


14 


0.031 


0.034 


15 


0.031 


0.034 


16 


0.031 


0.034 


17 


0.031 


0.033 


18 


0.031 


0.033 


19 


0.030 


0.033 


20 


0.031 


0.032 



Table 6.3 
Sphere Obstacle. 7 = 0.05. h = 1/32. 



Step 


\\E\\ L2 (n) 


B\\L 2 (n) 





0.383 


0.401 


1 


0.261 


0.264 


2 


0.176 


0.172 


3 


0.120 


0.113 


4 


0.081 


0.075 


5 


0.056 


0.051 


6 


0.039 


0.034 


7 


0.028 


0.025 


8 


0.021 


0.020 


9 


0.017 


0.017 


10 


0.015 


0.015 


11 


0.014 


0.015 


12 


0.013 


0.015 


13 


0.012 


0.015 


14 


0.013 


0.014 


15 


0.013 


0.014 


16 


0.013 


0.014 


17 


0.012 


0.015 


18 


0.013 


0.014 


19 


0.013 


0.014 


20 


0.013 


0.014 


Table 6.4 



Sphere Obstacle. 7 = 0.05. h = 1/64. 
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Fig. 6.1. Plot of Total Energy (\\E\\% 2 (fy + ll-^lli 2 (Q)/' vs - mes h s * ze after 20 time steps. The 
x-axis shows 1/h. 

7. Concluding remarks. We have shown that using appropriate finite-element 
spaces one can approximate the asymptotically-decaying solutions to Maxwell's equa- 
tions on the exterior of a spherical obstacle. We also have shown that the Crank- 
Nicholson time discretization keeps the electric and magnetic fields divergence free 
(the electric field weakly and the magnetic field-strongly). The next step is to apply 
the methods to more general domains and answer the question of whether ADS exist 
for obstacles with more complicated geometry. In relation to considering less regular 
obstacles, there are many open problems for the analysis of systems with dissipative 
boundary conditions. The computational techniques that we have introduced here 
easily generalize to cases of variable, matrix valued permittivity and permeability 
e(x) and n(x), as well as to more general hyperbolic systems for which the ADS 
phenomenon occurs. 
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